/*
** $Id: 3dmath.c,v 1.8 2006-01-14 05:07:12 xwyan Exp $
**
** Port to MiniGUI by Feynman from allegro.
**
** Copyright (C) 2006 Feynman Software.
 *
 * Vector and matrix manipulation routines.
 *
 * By Shawn Hargreaves.
 *
 * Allegro is a gift-ware.
**/

/*
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 2 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <minigui/common.h>
#include <minigui/minigui.h>
#include <minigui/gdi.h>
#include <minigui/window.h>
#include <minigui/control.h>
#include <minigui/fixedmath.h>

#include "mg3d.h"
#include "mg3dmatrix.h"
#include "mg3dmaths.h"


#ifndef MG_PI
   #define MG_PI   3.14159265358979323846
#endif

#define floatcos(x)     cos((x) * MG_PI / 128.0)
#define floatsin(x)     sin((x) * MG_PI / 128.0)
#define floattan(x)     tan((x) * MG_PI / 128.0)

mg3dMatrix identity_matrix = 
{
   {
      /* 3x3 identity */
      { 1<<16, 0,     0     },
      { 0,     1<<16, 0     },
      { 0,     0,     1<<16 },
   },

   /* zero translation */
   { 0, 0, 0 }
};

mg3dMatrixf identity_matrix_f = 
{
   {
      /* 3x3 identity */
      { 1.0, 0.0, 0.0 },
      { 0.0, 1.0, 0.0 },
      { 0.0, 0.0, 1.0 },
   },

   /* zero translation */
   { 0.0, 0.0, 0.0 }
};


/* scaling factors for the perspective projection */
fixed _persp_xscale = 160 << 16;
fixed _persp_yscale = 100 << 16;
fixed _persp_xoffset = 160 << 16;
fixed _persp_yoffset = 100 << 16;

float _persp_xscale_f = 160.0;
float _persp_yscale_f = 100.0;
float _persp_xoffset_f = 160.0;
float _persp_yoffset_f = 100.0;

#define CALC_ROW(n)     (fixmul(x, m->v[n][0]) +      \
	                 fixmul(y, m->v[n][1]) +      \
	                 fixmul(z, m->v[n][2]) +      \
	                        m->t[n])

void mg3dApplyMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z, fixed *xout, fixed *yout, fixed *zout)
{
    *xout = CALC_ROW(0);
    *yout = CALC_ROW(1);
    *zout = CALC_ROW(2);
}

#undef CALC_ROW

/* apply_matrix_f:
 *  Floating point vector by matrix multiplication routine.
 */
void mg3dApplyMatrixF(const mg3dMatrixf *m, float x, float y, float z, float *xout, float *yout, float *zout)
{
#define CALC_ROW(n) (x * m->v[(n)][0] + y * m->v[(n)][1] + z * m->v[(n)][2] + m->t[(n)])
       *xout = CALC_ROW(0);
       *yout = CALC_ROW(1);
       *zout = CALC_ROW(2);
#undef CALC_ROW
}


/*  mg3dGetTranslationMatrix:
 *  Constructs a 3d translation matrix. When applied to the vector 
 *  (vx, vy, vx), this will produce (vx+x, vy+y, vz+z).
 */
void mg3dGetTranslationMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z)
{
   *m = identity_matrix;

   m->t[0] = x;
   m->t[1] = y;
   m->t[2] = z;
}


fixed mg3dDotProduct(fixed x1, fixed y1, fixed z1, fixed x2, fixed y2, fixed z2)
{
   return fixmul(x1, x2) + fixmul(y1, y2) + fixmul(z1, z2);
}


float mg3dDotProductF(float x1, float y1, float z1, float x2, float y2, float z2)
{
   return (x1 * x2) + (y1 * y2) + (z1 * z2);
}


void mg3dPerspProject(fixed x, fixed y, fixed z, fixed *xout, fixed *yout)
{
   *xout = fixmul(fixdiv(x, z), _persp_xscale) + _persp_xoffset;
   *yout = fixmul(fixdiv(y, z), _persp_yscale) + _persp_yoffset;
}

void mg3dPerspProjectF (float x, float y, float z, float *xout, float *yout)
{
   float z1 = 1.0f / z;
   *xout = ((x * z1) * _persp_xscale_f) + _persp_xoffset_f;
   *yout = ((y * z1) * _persp_yscale_f) + _persp_yoffset_f;
}

/*  mg3dGetTranslationMatrixF:
 *  Floating point version of mg3dGetTranslationMatrix().
 */
void mg3dGetTranslationMatrixF(mg3dMatrixf *m, float x, float y, float z)
{
   *m = identity_matrix_f;

   m->t[0] = x;
   m->t[1] = y;
   m->t[2] = z;
}

/* mg3dGetScalingMatrix:
 *  Constructs a 3d scaling matrix. When applied to the vector 
 *  (vx, vy, vx), this will produce (vx*x, vy*y, vz*z).
 */
void mg3dGetScalingMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z)
{
   *m = identity_matrix;

   m->v[0][0] = x;
   m->v[1][1] = y;
   m->v[2][2] = z;
}

/*  mg3dGetScalingMatrixF:
 *  Floating point version of get_scaling_matrix().
 */
void mg3dGetScalingMatrixF(mg3dMatrixf *m, float x, float y, float z)
{
   *m = identity_matrix_f;

   m->v[0][0] = x;
   m->v[1][1] = y;
   m->v[2][2] = z;
}

/*  mg3dGetXRotMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around 
 *  the x axis by the specified amount (given in the Allegro fixed point, 
 *  256 degrees to a circle format).
 */
void mg3dGetXRotMatrix(mg3dMatrix *m, fixed r)
{
   fixed c = fixcos(r);
   fixed s = fixsin(r);

   *m = identity_matrix;

   m->v[1][1] = c;
   m->v[1][2] = -s;

   m->v[2][1] = s;
   m->v[2][2] = c;
}

/*  mg3dGetXRotMatrixF:
 *  Floating point version of get_x_rotate_matrix().
 */
void mg3dGetXRotMatrixF(mg3dMatrixf *m, float r)
{
   float c = floatcos(r);
   float s = floatsin(r);

   *m = identity_matrix_f;

   m->v[1][1] = c;
   m->v[1][2] = -s;

   m->v[2][1] = s;
   m->v[2][2] = c;
}

/*  mg3dGetYRotMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around 
 *  the y axis by the specified amount (given in the Allegro fixed point, 
 *  256 degrees to a circle format).
 */
void mg3dGetYRotMatrix(mg3dMatrix *m, fixed r)
{
   fixed c = fixcos(r);
   fixed s = fixsin(r);

   *m = identity_matrix;

   m->v[0][0] = c;
   m->v[0][2] = s;

   m->v[2][0] = -s;
   m->v[2][2] = c;
}

/*  mg3dGetYRotMatrixF:
 *  Floating point version of mg3dGetYRotMatrix().
 */
void  mg3dGetYRotMatrixF(mg3dMatrixf *m, float r)
{
   float c = floatcos(r);
   float s = floatsin(r);

   *m = identity_matrix_f;

   m->v[0][0] = c;
   m->v[0][2] = s;

   m->v[2][0] = -s;
   m->v[2][2] = c;
}

/*  mg3dGetZRotMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around 
 *  the z axis by the specified amount (given in the Allegro fixed point, 
 *  256 degrees to a circle format).
 */
void  mg3dGetZRotMatrix(mg3dMatrix *m, fixed r)
{
   fixed c = fixcos(r);
   fixed s = fixsin(r);

   *m = identity_matrix;

   m->v[0][0] = c;
   m->v[0][1] = -s;

   m->v[1][0] = s;
   m->v[1][1] = c;
}

/*  mg3dGetZRotMatrixF:
 *  Floating point version of mg3dGetZRotMatrix().
 */
void  mg3dGetZRotateMatrixF(mg3dMatrixf *m, float r)
{
   float c = floatcos(r);
   float s = floatsin(r);

   *m = identity_matrix_f;

   m->v[0][0] = c;
   m->v[0][1] = -s;

   m->v[1][0] = s;
   m->v[1][1] = c;
}

/* magical formulae for constructing rotation matrices */
#define MAKE_ROTATION(x, y, z)                  \
   fixed sin_x = fixsin(x);                       \
   fixed cos_x = fixcos(x);                       \
						\
   fixed sin_y = fixsin(y);                       \
   fixed cos_y = fixcos(y);                       \
						\
   fixed sin_z = fixsin(z);                       \
   fixed cos_z = fixcos(z);                       \
						\
   fixed sinx_siny = fixmul(sin_x, sin_y);        \
   fixed cosx_siny = fixmul(cos_x, sin_y);

#define MAKE_ROTATION_f(x, y, z)                \
   float sin_x = floatsin(x);                   \
   float cos_x = floatcos(x);                   \
						\
   float sin_y = floatsin(y);                   \
   float cos_y = floatcos(y);                   \
						\
   float sin_z = floatsin(z);                   \
   float cos_z = floatcos(z);                   \
						\
   float sinx_siny = sin_x * sin_y;             \
   float cosx_siny = cos_x * sin_y;

#define R00 (fixmul(cos_y, cos_z))
#define R10 (fixmul(sinx_siny, cos_z) - fixmul(cos_x, sin_z))
#define R20 (fixmul(cosx_siny, cos_z) + fixmul(sin_x, sin_z))

#define R01 (fixmul(cos_y, sin_z))
#define R11 (fixmul(sinx_siny, sin_z) + fixmul(cos_x, cos_z))
#define R21 (fixmul(cosx_siny, sin_z) - fixmul(sin_x, cos_z))

#define R02 (-sin_y)
#define R12 (fixmul(sin_x, cos_y))
#define R22 (fixmul(cos_x, cos_y))

#define R00_f (cos_y * cos_z)
#define R10_f ((sinx_siny * cos_z) - (cos_x * sin_z))
#define R20_f ((cosx_siny * cos_z) + (sin_x * sin_z))

#define R01_f (cos_y * sin_z)
#define R11_f ((sinx_siny * sin_z) + (cos_x * cos_z))
#define R21_f ((cosx_siny * sin_z) - (sin_x * cos_z))

#define R02_f (-sin_y)
#define R12_f (sin_x * cos_y)
#define R22_f (cos_x * cos_y)

/*  mg3dGetRotMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around
 *  all three axis by the specified amounts (given in the Allegro fixed 
 *  point, 256 degrees to a circle format).
 */
void  mg3dGetRotMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z)
{
   MAKE_ROTATION(x, y, z);

   m->v[0][0] = R00;
   m->v[0][1] = R01;
   m->v[0][2] = R02;

   m->v[1][0] = R10;
   m->v[1][1] = R11;
   m->v[1][2] = R12;

   m->v[2][0] = R20;
   m->v[2][1] = R21;
   m->v[2][2] = R22;

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dGetRotMatrixF:
 *  Floating point version of mg3dGetRotMatrix().
 */
void  mg3dGetRotMatrixF(mg3dMatrixf *m, float x, float y, float z)
{
   MAKE_ROTATION_f(x, y, z);

   m->v[0][0] = R00_f;
   m->v[0][1] = R01_f;
   m->v[0][2] = R02_f;

   m->v[1][0] = R10_f;
   m->v[1][1] = R11_f;
   m->v[1][2] = R12_f;

   m->v[2][0] = R20_f;
   m->v[2][1] = R21_f;
   m->v[2][2] = R22_f;

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dVectorLength: 
 *  Computes the length of a vector, using the son of the squaw...
 */
fixed  mg3dVectorLength(fixed x, fixed y, fixed z)
{
   x >>= 8;
   y >>= 8;
   z >>= 8;

   return (fixsqrt(fixmul(x,x) + fixmul(y,y) + fixmul(z,z)) << 8);
}

/*  mg3dVectorLengthF: 
 *  Floating point version of mg3dVectorLength().
 */
float  mg3dVectorLengthF(float x, float y, float z)
{
   return sqrt(x*x + y*y + z*z);
}

/*  mg3dNormalVector: 
 *  Converts the specified vector to a unit vector, which has the same
 *  orientation but a length of one.
 */
void  mg3dNormalVector(fixed *x, fixed *y, fixed *z)
{
   fixed length = mg3dVectorLength(*x, *y, *z);

   *x = fixdiv(*x, length);
   *y = fixdiv(*y, length);
   *z = fixdiv(*z, length);
}

/*  mg3dNormalVectorF: 
 *  Floating point version of  mg3dNormalVector().
 */
void  mg3dNormalVectorF(float *x, float *y, float *z)
{
   float length = 1.0 / mg3dVectorLengthF(*x, *y, *z);

   *x *= length;
   *y *= length;
   *z *= length;
}

/*  mg3dCrossProduct:
 *  Calculates the cross product of two vectors.
 */
void  mg3dCrossProduct(fixed x1, fixed y1, fixed z1, fixed x2, fixed y2, fixed z2, fixed *xout, fixed *yout, fixed *zout)
{
    *xout = fixmul(y1, z2) - fixmul(z1, y2);
    *yout = fixmul(z1, x2) - fixmul(x1, z2);
    *zout = fixmul(x1, y2) - fixmul(y1, x2);
}

/*  mg3dCrossProductF:
 *  Floating point version of mg3dCrossProduct().
 */
void  mg3dCrossProductF(float x1, float y1, float z1, float x2, float y2, float z2, float *xout, float *yout, float *zout)
{
    *xout = (y1 * z2) - (z1 * y2);
    *yout = (z1 * x2) - (x1 * z2);
    *zout = (x1 * y2) - (y1 * x2);
}

/*  mg3dGetAlignMatrix:
 *  Aligns a matrix along an arbitrary coordinate system.
 */
void  mg3dGetAlignMatrix(mg3dMatrix *m, fixed xfront, fixed yfront, fixed zfront, fixed xup, fixed yup, fixed zup)
{
   fixed xright, yright, zright;

    mg3dNormalVector(&xfront, &yfront, &zfront);
    mg3dNormalVector(&xup, &yup, &zup);

    mg3dCrossProduct(xfront, yfront, zfront, xup, yup, zup, &xright, &yright, &zright);
    mg3dCrossProduct(xright, yright, zright, xfront, yfront, zfront, &xup, &yup, &zup);

   m->v[0][0] = xright; 
   m->v[0][1] = xup; 
   m->v[0][2] = xfront; 

   m->v[1][0] = yright;
   m->v[1][1] = yup;
   m->v[1][2] = yfront;

   m->v[2][0] = zright;
   m->v[2][1] = zup;
   m->v[2][2] = zfront;

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dGetAlignMatrixF:
 *  Floating point version of mg3dGetAlignMatrix().
 */
void  mg3dGetAlignMatrixF(mg3dMatrixf *m, float xfront, float yfront, float zfront, float xup, float yup, float zup)
{
   float xright, yright, zright;

    mg3dNormalVectorF(&xfront, &yfront, &zfront);
    mg3dNormalVectorF(&xup, &yup, &zup);

    mg3dCrossProductF(xfront, yfront, zfront, xup, yup, zup, &xright, &yright, &zright);
    mg3dCrossProductF(xright, yright, zright, xfront, yfront, zfront, &xup, &yup, &zup);

   m->v[0][0] = xright; 
   m->v[0][1] = xup; 
   m->v[0][2] = xfront; 

   m->v[1][0] = yright;
   m->v[1][1] = yup;
   m->v[1][2] = yfront;

   m->v[2][0] = zright;
   m->v[2][1] = zup;
   m->v[2][2] = zfront;

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dGetVectorRotMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around
 *  the specified x,y,z vector by the specified angle (given in the Allegro 
 *  fixed point, 256 degrees to a circle format), in a clockwise direction.
 */
void  mg3dGetVectorRotMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z, fixed a)
{
   mg3dMatrixf rotation;
   int i, j;

   mg3dGetVectorRotMatrixF(&rotation, fixtof(x), fixtof(y), fixtof(z), fixtof(a));

   for (i=0; i<3; i++)
      for (j=0; j<3; j++)
	 m->v[i][j] = ftofix(rotation.v[i][j]);

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dGetVectorRotMatrixF:
 *  Floating point version of mg3dGetVectorRotMatrix().
 */
void  mg3dGetVectorRotMatrixF(mg3dMatrixf *m, float x, float y, float z, float a)
{
   float c = floatcos(a);
   float s = floatsin(a);
   float cc = 1 - c;

   mg3dNormalVectorF(&x, &y, &z);

   m->v[0][0] = (cc * x * x) + c;
   m->v[0][1] = (cc * x * y) + (z * s);
   m->v[0][2] = (cc * x * z) - (y * s);

   m->v[1][0] = (cc * x * y) - (z * s);
   m->v[1][1] = (cc * y * y) + c;
   m->v[1][2] = (cc * z * y) + (x * s);

   m->v[2][0] = (cc * x * z) + (y * s);
   m->v[2][1] = (cc * y * z) - (x * s);
   m->v[2][2] = (cc * z * z) + c;

   m->t[0] = m->t[1] = m->t[2] = 0;
}

/*  mg3dGetTransformMatrix:
 *  Constructs a 3d transformation matrix, which will rotate points around
 *  all three axis by the specified amounts (given in the Allegro fixed 
 *  point, 256 degrees to a circle format), scale the result by the
 *  specified amount (itofix(1) for no change of scale), and then translate
 *  to the requested x, y, z position.
 */
void  mg3dGetTransformMatrix(mg3dMatrix *m, fixed scale, fixed xrot, fixed yrot, fixed zrot, fixed x, fixed y, fixed z)
{
   MAKE_ROTATION(xrot, yrot, zrot);

   m->v[0][0] = fixmul(R00, scale);
   m->v[0][1] = fixmul(R01, scale);
   m->v[0][2] = fixmul(R02, scale);

   m->v[1][0] = fixmul(R10, scale);
   m->v[1][1] = fixmul(R11, scale);
   m->v[1][2] = fixmul(R12, scale);

   m->v[2][0] = fixmul(R20, scale);
   m->v[2][1] = fixmul(R21, scale);
   m->v[2][2] = fixmul(R22, scale);

   m->t[0] = x;
   m->t[1] = y;
   m->t[2] = z;
}

/*  mg3dGetTransformMatrixF:
 *  Floating point version of mg3dGetTransformMatrix().
 */
void  mg3dGetTransformMatrixF(mg3dMatrixf *m, float scale, float xrot, float yrot, float zrot, float x, float y, float z)
{
   MAKE_ROTATION_f(xrot, yrot, zrot);

   m->v[0][0] = R00_f * scale;
   m->v[0][1] = R01_f * scale;
   m->v[0][2] = R02_f * scale;

   m->v[1][0] = R10_f * scale;
   m->v[1][1] = R11_f * scale;
   m->v[1][2] = R12_f * scale;

   m->v[2][0] = R20_f * scale;
   m->v[2][1] = R21_f * scale;
   m->v[2][2] = R22_f * scale;

   m->t[0] = x;
   m->t[1] = y;
   m->t[2] = z;
}

/*  mg3dGetCameraMatrix: 
 *  Constructs a camera matrix for translating world-space objects into
 *  a normalised view space, ready for the perspective projection. The
 *  x, y, and z parameters specify the camera position, xfront, yfront,
 *  and zfront is an 'in front' vector specifying which way the camera
 *  is facing (this can be any length: normalisation is not required),
 *  and xup, yup, and zup is the 'up' direction vector. Up is really only
 *  a 1.5d vector, since the front vector only leaves one degree of freedom
 *  for which way up to put the image, but it is simplest to specify it
 *  as a full 3d direction even though a lot of the information in it is
 *  discarded. The fov parameter specifies the field of view (ie. width
 *  of the camera focus) in fixed point, 256 degrees to the circle format.
 *  For typical projections, a field of view in the region 32-48 will work
 *  well. Finally, the aspect ratio is used to scale the Y dimensions of
 *  the image relative to the X axis, so you can use it to correct for
 *  the proportions of the output image (set it to 1 for no scaling).
 */
void  mg3dGetCameraMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z, fixed xfront, fixed yfront, fixed zfront, fixed xup, fixed yup, fixed zup, fixed fov, fixed aspect)
{
   mg3dMatrixf camera;
   int i, j;

   mg3dGetCameraMatrixF(&camera,
		       fixtof(x), fixtof(y), fixtof(z), 
		       fixtof(xfront), fixtof(yfront), fixtof(zfront), 
		       fixtof(xup), fixtof(yup), fixtof(zup), 
		       fixtof(fov), fixtof(aspect));

   for (i=0; i<3; i++) {
      for (j=0; j<3; j++)
	 m->v[i][j] = ftofix(camera.v[i][j]);

      m->t[i] = ftofix(camera.t[i]);
   }
}

/*  mg3dGetCameraMatrixF: 
 *  Floating point version of mg3dGetCameraMatrix().
 */
void  mg3dGetCameraMatrixF(mg3dMatrixf *m, float x, float y, float z, float xfront, float yfront, float zfront, float xup, float yup, float zup, float fov, float aspect)
{
   mg3dMatrixf camera, scale;
   float xside, yside, zside, width, d;

   /* make 'in-front' into a unit vector, and negate it */
   mg3dNormalVectorF(&xfront, &yfront, &zfront);
   xfront = -xfront;
   yfront = -yfront;
   zfront = -zfront;

   /* make sure 'up' is at right angles to 'in-front', and normalize */
   d =  mg3dDotProductF(xup, yup, zup, xfront, yfront, zfront);
   xup -= d * xfront; 
   yup -= d * yfront; 
   zup -= d * zfront;
   mg3dNormalVectorF(&xup, &yup, &zup);

   /* calculate the 'sideways' vector */
    mg3dCrossProductF(xup, yup, zup, xfront, yfront, zfront, &xside, &yside, &zside);

   /* set matrix rotation parameters */
   camera.v[0][0] = xside; 
   camera.v[0][1] = yside;
   camera.v[0][2] = zside;

   camera.v[1][0] = xup; 
   camera.v[1][1] = yup;
   camera.v[1][2] = zup;

   camera.v[2][0] = xfront; 
   camera.v[2][1] = yfront;
   camera.v[2][2] = zfront;

   /* set matrix translation parameters */
   camera.t[0] = -(x*xside  + y*yside  + z*zside);
   camera.t[1] = -(x*xup    + y*yup    + z*zup);
   camera.t[2] = -(x*xfront + y*yfront + z*zfront);

   /* construct a scaling matrix to deal with aspect ratio and FOV */
   width = floattan(64.0 - fov/2);
    mg3dGetScalingMatrixF(&scale, width, -aspect*width*4/3, -1.0);

   /* combine the camera and scaling matrices */
    mg3dMatrixMulF(&camera, &scale, m);
}

/*  mg3dQTranslateMatrix:
 *  Adds a position offset to an existing matrix.
 */
void  mg3dQTranslateMatrix(mg3dMatrix *m, fixed x, fixed y, fixed z)
{
   m->t[0] += x;
   m->t[1] += y;
   m->t[2] += z;
}

/*  mg3dQTranslateMatrixF:
 *  Floating point version of mg3dQTranslateMatrix().
 */
void  mg3dQTranslateMatrixF(mg3dMatrixf *m, float x, float y, float z)
{
   m->t[0] += x;
   m->t[1] += y;
   m->t[2] += z;
}

/*  mg3dQscaleMatrix:
 *  Adds a scaling factor to an existing matrix.
 */
void  mg3dQscaleMatrix(mg3dMatrix *m, fixed scale)
{
   int i, j;

   for (i=0; i<3; i++)
      for (j=0; j<3; j++)
	 m->v[i][j] = fixmul(m->v[i][j], scale);
}

/*  mg3dQscaleMatrixF:
 *  Floating point version of mg3dQscaleMatrix().
 */
void  mg3dQscaleMatrixF(mg3dMatrixf *m, float scale)
{
   int i, j;

   for (i=0; i<3; i++)
      for (j=0; j<3; j++)
	 m->v[i][j] *= scale;
}

/*  mg3dMatrixMul:
 *  Multiplies two matrices, storing the result in out (this must be
 *  different from the two input matrices). The resulting matrix will
 *  have the same effect as the combination of m1 and m2, ie. when
 *  applied to a vector v, (v * out) = ((v * m1) * m2). Any number of
 *  transformations can be concatenated in this way.
 */
void  mg3dMatrixMul(const mg3dMatrix *m1, const mg3dMatrix *m2, mg3dMatrix *out)
{
   mg3dMatrix temp;
   int i, j;

   if (m1 == out) {
      temp = *m1;
      m1 = &temp;
   }
   else if (m2 == out) {
      temp = *m2;
      m2 = &temp;
   }

   for (i=0; i<3; i++) {
      for (j=0; j<3; j++) {
	 out->v[i][j] = fixmul(m1->v[0][j], m2->v[i][0]) +
			fixmul(m1->v[1][j], m2->v[i][1]) +
			fixmul(m1->v[2][j], m2->v[i][2]);
      }

      out->t[i] = fixmul(m1->t[0], m2->v[i][0]) +
		  fixmul(m1->t[1], m2->v[i][1]) +
		  fixmul(m1->t[2], m2->v[i][2]) +
		  m2->t[i];
   } 
}

/*  mg3dMatrixMulF:
 *  Floating point version of mg3dMatrixMul().
 */
void  mg3dMatrixMulF(const mg3dMatrixf *m1, const mg3dMatrixf *m2, mg3dMatrixf *out)
{
   mg3dMatrixf temp;
   int i, j;

   if (m1 == out) {
      temp = *m1;
      m1 = &temp;
   }
   else if (m2 == out) {
      temp = *m2;
      m2 = &temp;
   }

   for (i=0; i<3; i++) {
      for (j=0; j<3; j++) {
	 out->v[i][j] = (m1->v[0][j] * m2->v[i][0]) +
			(m1->v[1][j] * m2->v[i][1]) +
			(m1->v[2][j] * m2->v[i][2]);
      }

      out->t[i] = (m1->t[0] * m2->v[i][0]) +
		  (m1->t[1] * m2->v[i][1]) +
		  (m1->t[2] * m2->v[i][2]) +
		  m2->t[i];
   }
}

/*  mg3dPolygonZNormal:
 *  Helper function for backface culling: returns the z component of the
 *  normal vector to the polygon formed from the three vertices.
 */
fixed  mg3dPolygonZNormal(const mg3dVp *v1, const mg3dVp *v2, const mg3dVp *v3)
{
   return (fixmul(v2->x-v1->x, v3->y-v2->y) - fixmul(v3->x-v2->x, v2->y-v1->y));
}

/*  mg3dPolygonZNormalF:
 *  Floating point version of mg3dPolygonZNormal().
 */
float  mg3dPolygonZNormalF(const mg3dVpf *v1, const mg3dVpf *v2, const mg3dVpf *v3)
{
   return ((v2->x-v1->x) * (v3->y-v2->y)) - ((v3->x-v2->x) * (v2->y-v1->y));
}

/*  mg3dSetProjectionViewport:
 *  Sets the viewport used to scale the output of the mg3dPerspProject() 
 *  function.
 */
void  mg3dSetProjectionViewport(int x, int y, int w, int h)
{
   _persp_xscale = itofix(w/2);
   _persp_yscale = itofix(h/2);
   _persp_xoffset = itofix(x + w/2);
   _persp_yoffset = itofix(y + h/2);

   _persp_xscale_f = w/2;
   _persp_yscale_f = h/2;
   _persp_xoffset_f = x + w/2;
   _persp_yoffset_f = y + h/2;
}


